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Abstract  —  In  reservoir  simulation,  the  most  known  well-reservoir  coupling  technique  is  based  on  the 
Peaceman’s  equivalent  radius,  which  is  applied  on  the  productivity  index  computation,  that  is  used  to 
relate  flow  rate,  wellbore  pressure,  and  the  mesh  block  pressure.  Original  Peaceman’s  model  considers 
the  assumption  of  steady-state  flow,  leading  to  an  artifact  on  the  calculation  of  the  wellbore  pressure, 
called  numerical  wellbore  storage.  This  artifact  is  more  significant  for  coarser  meshes  and  initial  time 
instants,  and  it  is  also  a  function  of  fluid  and  rock  properties.  In  this  context,  we  have  implemented  and 
compared  different  Peaceman’s  technique  extensions  for  productivity  index  calculation  incorporating 
transient  effects  to  prevent  numerical  wellbore  storage.  We  have  also  considered  some  production 
scenarios  for  a  vertical  well  and  single-phase  oil  flow. 

Keywords — Numerical  storage,  Peaceman’s  equivalent  radius,  Transient  well  index,  Well-reservoir 
coupling. 


I.  INTRODUCTION 


of  the  fluid  produced  at  the  surface.  Knowledge  of  the 
behavior,  under  pressure  and  temperature  variations, 
of  oii,  naturai  gas,  and  water,  aione  or  in  combina¬ 
tion,  is  of  fundamentai  interest  to  Petroieum  Engineer¬ 
ing  [14],  whether  under  static  or  moving  conditions,  in 
the  rocks  of  the  reservoir  or  ducts. 


During  the  1960s,  the  term  Numericai  Simuiation 
of  Reservoirs  became  common  in  the  oii  industry  [14], 
referring  to  the  use  of  physicai-mathematicai  modeis 
and  computationai  toois  to  predict  the  performance 
of  hydrocarbon  reservoirs  under  different  production 
scenarios.  According  to  Dumkwu  et  ai.  [16],  Numeri¬ 
cai  Reservoir  Simuiation  has  become  an  aiiy  in  Reser¬ 
voir  Engineering,  and  we  use  it  in  decision  making  that 
invoives  many  financiai  resources,  in  estimating  un¬ 
derground  reserves,  and  to  forecast  and  improve  per¬ 
formance  production  of  reservoirs  in  the  oii  and  gas 
industry.  This  work  focuses  on  the  caicuiation  of  pres¬ 
sure  estimates  in  weiis,  using  weii-reservoir  coupiing 
techniques. 


Over  the  years,  Petroieum  Engineering  has  rec¬ 
ognized  the  need  for  accurate  information  reiated  to 
physicai  conditions  in  the  weii  and  inside  the  reser¬ 
voir.  With  the  initiai  progress  in  oii  recovery  methods, 
it  became  ciear  that  caicuiations  made  with  informa¬ 
tion  from  the  surface  or  the  top  of  the  weii  couid  of¬ 
ten  iead  to  misunderstandings.  Sciater  et  ai.  [33]  de¬ 
scribed  the  first  pressure  vaiue  record  at  the  bottom 
of  the  weii  and  fluid  coiiection  in  the  weii,  under  pres¬ 
sure,  for  sampiing.  These  same  authors  define  rock 
bottom  data  by  referring  to  pressure,  temperature, 
gas-oii  ratio,  and  the  chemicai  and  physicai  nature 
of  the  fluid.  Miiiikan  et  ai.  [21]  highiighted  the  need 
of  accurate  weiibore  pressure  measurements  when 
they  described  the  first  accurate  pressure  gauge,  and 
they  pointed  to  the  fundamentai  importance  of  weii¬ 
bore  pressure  knowiedge  for  Petroieum  Engineering, 


1 .1 .  Determination  of  pressure  in  oii  weiis 


Petroieum  is  the  name  given  to  naturai  hydrocar¬ 
bon  mixtures  that  can  be  found  in  the  soiid,  iiquid,  and 
gas  phases,  depending  on  the  pressure  and  temper¬ 
ature  conditions.  In  nature,  oil  can  appear  in  only  one 
phase,  or  it  can  appear  as  a  multiphase  mixture  in 
equilibrium  [30].  In  many  cases,  the  fluid  phase  in  the 
reservoir  depends  on  the  thermodynamics  properties 
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in  search  of  more  efficient  recovery  and  iifting  proce¬ 
dures.  Based  on  this  contribution,  we  were  abie  to 
measure  pressure,  which  is  essentiai  information  for 
reservoir  performance  caicuiations  [14], 

The  Weii  Testing  Anaiysis  is  the  engineering 
branch  which  obtains  estimates  of  properties  of  the 
weii-reservoir  system  from  weiibore  data  under  pro¬ 
duction/injection  conditions,  as  weii  as  through  the 
appiication  of  inverse  probiems.  We  are  interested 
in  the  pressure  variation  in  the  weii  as  a  function  of 
time.  We  accompiish  this  by  measuring  the  pressures 
at  the  bottom  of  the  weii  and  the  flow  at  the  surface, 
for  exampie.  From  the  measured  pressure  response, 
it  is  possibie  to  determine  reservoir  properties  usefui 
for  pianning  the  compietion  of  the  weii  or  the  reser¬ 
voir  depietion  pian  [35].  During  a  weii  test,  we  can 
create  a  transient  pressure  response  by  a  temporary 
change  in  the  production  rate.  The  weii  response  is 
usuaiiy  monitored  for  a  reiativeiy  short  period,  com¬ 
pared  to  the  reservoir  iife,  depending  on  the  test  ob¬ 
jectives.  For  estimates  invoiving  regions  around  the 
weiis,  we  often  carried  out  tests  in  iess  than  two  days. 
In  the  case  of  tests  to  analyze  the  reservoir  bound¬ 
aries,  it  may  take  several  months  to  register  pressure 
data  [20]. 

It  is  usual  to  use  analytical  solutions  to  interpret 
the  results  of  pressure  tests  in  wells.  Theis  [37]  pro¬ 
posed  one  of  the  first  theoretical  solutions  for  deter¬ 
mining  the  well  pressure,  depending  on  the  flow  rate 
and  the  duration  of  production.  Considering  the  one¬ 
dimensional  flow  in  the  radial  direction  in  a  reservoir 
described  in  Cylindrical  coordinates,  producing  at  a 
constant  flow  Qsc  (defined  in  standard  conditions,  neg¬ 
ative  for  production  and  positive  for  injection),  from  the 
initial  time  t=0,  the  analytical  solution  for  the  pressure 
in  the  reservoir  in  a  radius  r  is  given  by  [30]: 


p{r,  t)  =pi  + 


qscBp 

Airkh 


4>iiCtr‘^  \ 
Akt  j\  ’ 


(1) 


where  pi  is  the  initial  pressure  of  the  porous  medium, 
fj.  is  the  viscosity  of  the  fluid,  k  is  the  permeability 
of  the  medium,  h  is  the  reservoir  thickness,  Ei  is 
the  Integral  Exponential  Function,  4>  is  the  porosity  of 
the  medium  and  ct  is  the  total  compressibility  that  in¬ 
cludes  the  contribution  of  the  fluid  (co)  and  rock  (c^) 
compressibilities  (c*  =  Co  +  c^). 

For  small  values  of  the  argument  (X  <  0.025), 
the  integral  exponential  function  can  be  approximated. 


with  an  error  of  less  than  1%,  by  [26], 


E,{-X)  ^  ln{XX), 


(2) 


where  A  =  exp(0. 57722),  and  0.57722  is  the  Euler 
constant.  Thus,  the  pressure  in  the  well  (p^,/)  is  ap¬ 
proximated  by 


Pwf{t)  =Pi  + 


qscBp 

Ankh 


f  X(j)pctrl,y 

V  j\  ’ 


(3) 


where  is  the  radius  of  the  well. 

Throughout  the  years,  researchers  have  devel¬ 
oped  models  for  more  complex  reservoirs.  However, 
there  are  many  situations  for  which  analytical  solu¬ 
tions  do  not  exist  or  are  difficult  to  obtain.  Hetero¬ 
geneous  reservoirs  with  irregular  borders  and  with 
complex  well  geometry  are  examples  where  numeri¬ 
cal  simulation  is  the  alternative  to  determine  transient 
pressure  responses. 


1.2.  Numerical  reservoir  simulation 

As  already  said,  we  often  use  reservoir  simulation 
models  in  the  oil  and  gas  industry.  The  reasons  for 
such  acceptance  are  due  to  some  advances  [17]: 

1 .  in  computing  (particularly  concerning  the  speed 
of  processing  and  the  increase  in  the  storage 
capacity  of  computers); 

2.  numerical  techniques  for  solving  partial  differen¬ 
tial  equations  (PDEs); 

3.  reservoir  simulators  suitable  for  use  in  modeling 
field  cases; 

4.  reservoir  characterization  techniques,  among 
other  factors. 

The  ultimate  goal  of  a  numerical  simulation  study 
is  to  accurately  predict  flow  and  pressure  in  the  well 
and  estimate  pressure  distributions  in  the  reservoir 
(and  saturation,  in  the  multiphase  case).  When  we  in¬ 
tend  to  incorporate  well  modeling  in  reservoir  simula¬ 
tion,  some  difficulties  appear,  which  we  can  separate 
into  three  main  categories  [17]: 

1 .  The  cell  of  the  computational  mesh  that  contains 
the  well  is  usually  large  compared  to  the  well¬ 
bore  radius,  due  to  the  discrepancy  between  the 
spatial  scales  involved.  The  average  diameter 
of  the  production/injection  well  is  about  10  cm, 
while  the  dimensions  of  a  reservoir  can  be  kilo¬ 
meters.  Therefore,  the  cell  pressure  calculated 
is  a  poor  estimate  of  the  pressure  in  the  well; 
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2.  Coupling  is  a  complex  interaction  between  the 
reservoir  (porous  medium)  and  the  well  (a 
duct),  particularly  in  the  case  of  wells  that  pass 
through  several  layers  of  the  porous  medium; 

3.  There  are  difficulties  in  describing  the  well- 
reservoir  coupling  in  the  multiphase  flow  when 
the  total  production  flow  of  the  well  is  specified. 

It  is  also  worth  noting  that  there  is  a  variation  in 
the  time  scale,  with  phenomena  involving  strong  pres¬ 
sure  gradients  that  may  occur  in  short  periods,  such 
as  fractions  of  an  hour,  and  border  effects  (such  as,  for 
example,  pressure  drop  at  the  reservoir  boundaries 
concerning  the  initial  pressure)  that  may  take  months 
to  appear.  These  times  are  a  function  of  the  charac¬ 
teristics  of  the  well-reservoir  system.  Other  problems 
can  also  arise  when  considering  the  location  of  the 
producer/injector  well  displaced  from  the  center  of  the 
cell,  or  even  when  there  are  several  wells  in  a  cell  in 
the  computational  mesh. 

The  well-reservoir  coupling  model  appears  as  a 
central  component  in  reservoir  simulation.  Peaceman 
[26]  introduced  a  well-known  and  commonly  used 
model,  which  considers  that  the  flow  between  the 
reservoir  and  the  producing  well  is  steady-state.  Ex¬ 
tensions  were  also  created,  for  example,  for  the  case 
of  anisotropy  [27],  the  incorporation  of  transient  ef¬ 
fects  [9],  and  techniques  for  horizontal  wells  [2].  Al¬ 
though  this  coupling  model  provides  suitable  results 
for  reservoirs  in  a  wide  range  of  applications,  including 
their  extensions,  it  introduces  a  phenomenon  called 
numerical  storage  (due  to  the  similarity  of  the  results 
when  considering  wellbore  storage,  although  it  does 
not  correspond  to  the  physical  behavior  in  the  well) 
for  the  initial  instants,  and  the  application  of  method¬ 
ologies  to  circumvent  this  numerical  issue  is  the  main 
objective  of  this  work. 

Therefore,  this  work  aims  to  implement  a  model 
for  a  transient  coupling  of  the  well-reservoir  system  to 
minimize  the  phenomenon  of  numerical  storage  orig¬ 
inating  from  the  steady-state  well-reservoir  coupling, 
using  the  technique  of  Peaceman  [26].  We  consider 
here  a  single-phase  flow,  of  a  slightly  compressible  oil, 
in  a  reservoir  containing  a  vertical  producer  well  that 
penetrates  the  entire  oil  reservoir.  We  also  use  a  poly¬ 
nomial  function,  with  low  error  and  wide  application 
range,  when  calculating  the  integral  exponential  Et. 


II.  FLOW  MODELING  IN  POROUS  MEDIA 

In  this  work,  we  consider  the  following  hypotheses 
for  the  flow  in  the  porous  medium: 

1.  the  porous  medium  is  homogeneous  and 
anisotropic  in  terms  of  its  permeability; 

2.  porosity  is  a  linear  function  of  pressure; 

3.  the  compressibility  of  the  rock  is  small  and  con¬ 
stant; 

4.  the  fluid  is  Newtonian  and  slightly  compressible; 

5.  the  fluid  has  a  constant  chemical  composition; 

6.  there  are  no  electrokinetic  effects; 

7.  there  are  no  inertial  or  turbulent  effects; 

8.  the  flow  is  single-phase  and  isothermal; 

9.  there  are  no  chemical  reactions. 

2.1.  Governing  equations 

The  continuity  equation  is  a  mass  balance  equa¬ 
tion  so  that  the  difference  in  mass  entering  and  leav¬ 
ing  a  control  volume  must  be  equal  to  the  accumula¬ 
tion  of  mass  within  the  control  volume  in  a  given  time 
interval.  In  oil  reservoirs,  the  control  volume  is  a  por¬ 
tion  of  the  porous  medium  that  can  contain  one,  two, 
or  three  phases  of  fluid  [17].  It  is  worth  mentioning 
that  we  consider  the  porous  medium  as  a  continuum. 
So,  we  define  the  physical  properties,  at  any  point,  us¬ 
ing  the  concept  of  Representative  Elementary  Volume 
(REV)  [41].  For  single-phase  mass  flows  in  a  porous 
medium,  we  can  write  the  continuity  equation  in  the 
form: 

^  +  V  •  (pv)  -  ^  =  0  (4) 

where  p  is  the  density  of  the  fluid,  v  is  the  apparent 
fluid  velocity  (flow  rate  divided  by  the  cross-sectional 
area),  qm  is  a  source  term  representing  the  production 
or  injection  of  fluid,  and  14  is  the  volume  of  the  rock 
(solid  material  plus  pores). 

We  can  also  write  the  continuity  equation  in  terms 
of  the  formation  volume  factor  [17],  B  =  Psc/p,  as 


where  the  sc  subscript  indicates  the  standard  condi¬ 
tions. 

The  civil  engineer  Henry  Darcy,  in  1856,  through 
experiments  of  vertical  water  filtration  in  columns  of 
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homogeneous  sand,  obtained  the  first  equation  used 
to  describe  the  movement  of  fluids  in  porous  media. 
This  equation  was  iater  caiied  Darcy’s  iaw.  It  is  an 
empiricai  reiationship  between  the  mass  flow  rate  of 
fluid  through  a  porous  medium  and  the  potentiai  gra¬ 
dient,  defined  in  the  form  [15] 

-V$  =  ^,  (6) 

where  K  symboiizes  the  hydrauiic  conductivity  of  the 
porous  medium,  which  combines  fluid  and  medium 
properties,  and  V<i>  is  the  potentiai  gradient,  defined 
as  [17,  30], 

V$  =  Vp  — 7gV£)  (7) 


where  jg  =  pg,  g  is  the  magnitude  of  the  acceieration 
due  to  gravity,  and  D  is  the  depth  (positive  in  the  ver- 
ticaiiy  downward  direction).  Later,  Darcy’s  iaw  couid 
be  deduced  mathematicaiiy  from  the  baiance  equa¬ 
tions  that  govern  the  singie-phase  flow  of  a  Newtonian 
fluid  [40]. 

Rewriting  Darcy’s  iaw  in  terms  of  the  properties  of 
the  fluid,  we  can  express  the  surface  flow  velocity  as 


V  =  --  {Vp-  pgVD) , 
P 


(8) 


where  k  is  the  absolute  permeability  tensor. 

From  the  definition  of  the  compressibility  coeffi¬ 
cients  of  the  rock  and  the  fluid,  it  is  possible  to  rewrite 
the  transient  term  of  Eq.  (5)  in  the  form 


dt\B 


d  /  1\  dp  1  d(j)  dp 
dp\B)  dt^Bdpdt 


i- 

dp 


1  d(j) 
B  dp 


dp 


4>Co 

B° 


B 


dp 

dt' 


(9) 


where  we  considered  that  (p  =  (j){p)  and  B  =  B{p)  and 
that  [17] 


0  =  [1  +  C0(p-p°)]  , 


(10) 


P  =  [l  +  Co(p-P°)]  , 


(11) 


B° 

[1  +  Co{p  -  p°)]’ 


(12) 


ri  r  OM  ’  (1  ^) 

[1  -  Cf,{p-p^)] 

where  the  superscript  0  indicates  the  reference  con¬ 
ditions,  and  is  the  coefficient  of  variation  of  viscos¬ 
ity,  and  we  assume  that  Co  and  are  small  and  con¬ 
stants. 

We  obtain  the  Hydraulic  Diffusivity  Equation  (HDE) 
by  combining  Darcy’s  law  with  the  mass  conservation 
equation.  Now,  writing  Eq.  (8)  in  terms  of  its  compo¬ 
nents. 


Vx  = 


p 


p 


Vz  =  —  \  - 
P 


dp  dD\ 

)  ’ 

(14) 

dp  dD\ 

dy'^  dy  )  ' 

(15) 

dp  dD\ 

(16) 

In  turn,  we  can  write  the  conservation  equation  as 


d 

di 


B 


JL 

dx 


d  /Vy\  a 

Vs/  dyKBJ  dz\B) 


dz 


Qm 

^bPsc 

(17) 


Next,  we  multiply  Eq.  (17)  by  V),  =  dxdydz, 


d  f  (p 


dt\B 


d_ 

dx 

_9 

dz 


'^x^x 

B 

VzAz 

B 


dx  — 


d_ 

dy 


'^yAy 

B 


j  I 

az  H - , 


dy 


(18) 


where  Ax,  Ay,  and  A^  are,  respectively,  the  areas  of 
surfaces  normal  to  the  x-,  y-  and  z-  directions. 

Then,  replacing  Vx,  Vy,  and  Vz  and  using  the  defi¬ 
nition  Qsc  =  Qm/Psc  we  finally  obtain 


d  ( (p 


^'’dt  I 


A 

dx 
d 

+  TT- 


Ax  kx 
pB 
A  k 

J-lylx,y 


dy  [  pB 


-F 


-F 


d 


Azkz 


dz  [  pB 

Qsc- 


dp 

dx 

dp 

dy 

dp 

dz 


dD 


pg 


pg 


pg 


dx 

dy  ) 

dD 

dz 


dx 


dy 


dz 


(19) 


Now,  using  the  result  obtained  in  Eq.  (9), 


dp 


d_ 

dx 

d_ 

dy 

dz 


Axkx 
pB 
A  h 
pB 

Azkz 

pB 


dp 

dx 

dp 

dy 

dp 

dz 


pg 


pg 


pg 


dD 

dx 

dD 

dy 

dD 

dz 


dx 

dy 

dz  -F  qac 


(20) 
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where  the  Fp  coefficient  represents  the  compressibii- 
ity  effects  of  rock  and  fluid,  given  by 


rp  =  Vb 


B  ) 


(21) 


Incorporating  the  gravitationai  effects  in  the  term 
Fg,  we  can  rewrite  the  HDE  as  foiiows: 

dp  d  ( dp\  ,  d  f  Ayky  dp\ 

d  /  A^kz  dp\  ,  ^ 

where 


Tg  = 


+ 


d^ 

dx 

d_ 

Wz 


-P9-^  I  dx 

OX 

I  dy. 


dy 


(23) 


2.2.  Initiai  and  boundary  conditions 

To  soive  Eq.  (22),  we  need  to  provide  the  appropri¬ 
ate  initiai  and  boundary  conditions.  As  an  initiai  con¬ 
dition,  we  take 

p{x,y,z,t  =  0)=p,niix,y,z),  (24) 

where  k™  is  the  initiai  pressure  distribution  before  the 
reservoir  is  disturbed  by  fluid  production/injection. 

As  externai  boundary  conditions  (externai  borders 
of  the  reservoir),  we  can  use  prescribed  pressure  or 
flow  rate.  For  a  sealed  rectangular  parallelepiped 
reservoir  of  edges  Ly  and  L^,  we  have 

^  ^  ^ 
\C>y)y=o,Ly  U^/z=0.L. 

(25) 


how  we  can  use  expressions  of  this  type,  we  assume 
steady-state  flow  in  the  region  close  to  the  well.  Fur¬ 
ther,  we  consider  the  radial  flow  of  an  incompressible 
fluid  towards  the  vertical  well  of  radius  r^,,  in  a  rock  for¬ 
mation  with  uniform  permeability  and  thickness.  So, 
assuming  these  conditions,  it  is  possible  to  write  [17]: 


—2nkHhr  dp 
^  pB  Wr 


(26) 


where  kn  represents  the  permeability  value  in  the  ra¬ 
dial  direction. 

Integrating  Eq.  (26)  between  the  radius  of  the  well, 
Tu,,  and  an  arbitrary  radius,  r,  where  <  r  <  re  (re 
is  the  external  radius),  it  is  possible  to  obtain  an  equa¬ 
tion  for  the  pressure  distribution  for  the  steady-state 
well-coupling. 


P=Pwf  - 


2TTkHh 


(27) 


For  r  =  re,  where  we  define  the  pressure  as  pe,  we 
can  rewrite  Eq.  (27)  as 


—2'KkHh 
pBln  (  — 


Pwf)  : 


(28) 


so  that  Eq.  (28)  provides  the  well  production  in  terms 
of  external  and  well  pressures. 

In  general,  we  express  the  production  rate,  under 
standard  conditions,  using  the  pressure  of  the  well 
and  the  average  pressure  of  the  reservoir  [17].  From 
the  development  of  Van  et  al.  [39],  for  a  reservoir  de¬ 
scribed  in  Cylindrical  coordinates,  considering  the  av¬ 
erage  volumetric  pressure  in  the  oil  reservoir,  p,  be¬ 
tween  Tyj  and  Te,  and  even  though  re  >  Tw,  we  have 


2.3.  Source  term 

In  the  numerical  simulation  of  reservoirs,  wells  are 
considered  internal  boundaries,  and,  in  the  case  of 
the  use  of  Cartesian  coordinates,  the  term  source, 
used  to  represent  wells,  is  directly  related  to  the  use 
of  well-reservoir  coupling  techniques.  We  can  impose 
these  boundary  conditions  specifying  the  pressure  in 
the  well  (Dirichlet  condition)  or  the  flow  rate  in  the  well 
(Neumann  condition)  [17]. 

For  the  representation  of  wells,  it  is  necessary  to 
obtain  an  expression  that  relates  the  pressure  in  the 
porous  medium,  p,  with  the  pressure  in  the  well,  p^f, 
and  the  flow  rate  in  the  well,  As  an  example  of 


Qsc  — 


—2'KkHh 


pB 


(p-Pwf)  , 


(29) 


which  differs  from  Eq.  (28)  only  by  the  skin  factor  s, 
the  1/2  in  the  denominator,  and  by  replacing  pe  with  p. 
The  introduction  of  the  skin  factor  makes  Eq.  (29) 
more  general,  valid  for  the  case  where  there  is  an  ad¬ 
ditional  pressure  drop  around  the  well,  for  example, 
due  to  the  formation  damage.  For  the  specific  case 
in  which  the  permeability  may  change,  we  determine 
the  s  factor  by  [18,  30] 


(30) 
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where  ks  is  the  modified  permeabiiity,  and  rs  is  the 
radius  of  the  region  where  the  permeabiiity  can  vary. 
When  ks  <  k,  the  vaiue  of  s  is  positive,  and  we  have 
damage  to  the  formation.  Otherwise,  s  is  negative, 
and  we  have  stimuiation.  On  the  other  hand,  when 
s  =  0,  there  is  no  influence  on  the  productivity  of  the 
weii  by  iocai  permeabiiity  changes  near  the  weii  [30]. 

Commonly,  we  rewrite  Eq.  (29)  in  a  more  compact 
form,  introduced  by  Van  et  al.  [39], 


Qsc  -  Jwi/P  Pwf  ) 


where  is  the  Productivity  Index  (PI): 


2'KkHh 


fj.B 


+  s  —  F 


(31) 


(32) 


or  yet, 

^  (33) 

where  is  the  geometric  factor  of  the  well,  defined 
as 


Gu, 


2Trkh 


+  s  —  F 


(34) 


The  factor  F,  in  Eqs.  (32)  and  (34),  assumes  the 
values  of  1/2  for  steady-state  conditions  and  3/4  for 
the  pseudo-steady  state  [17].  The  representation  of 
Eq.  (32)  for  the  Productivity  Index  comes  from  ana¬ 
lytical  solutions.  However,  we  also  apply  the  concept 
of  the  Productivity  Index  in  the  context  of  numerical 
reservoir  simulation,  and  its  calculation  is  associated, 
in  this  case,  with  the  use  of  well-reservoir  coupling 
techniques. 


III.  NUMERICAL  METHODOLOGY 

We  describe  fluid  flow  in  porous  media  using  a  set 
of  partial  differential  equations,  and  we  can  not  always 
obtain  analytical  solutions  due  to  the  non-linear  nature 
of  the  equations.  Therefore,  we  must  use  numerical 
techniques  to  solve  the  balance  equations  that  govern 
the  flow.  Among  the  numerical  methods  that  we  can 
apply  in  its  resolution,  we  choose  to  use  the  Finite  Dif¬ 
ference  Method  (EDM),  which  is  the  most  used  in  the 
oil  industry  [17,  26]. 


3.1.  Discretization 

Discretization  is  the  process  of  converting  the 
PDE,  valid  in  the  continuous  medium  that  defines  the 
resolution  domain,  by  a  set  of  algebraic  equations  de¬ 
fined  in  a  discrete  domain,  to  obtain  a  numerical  so¬ 
lution.  The  first  step  in  the  discretization  stage  is  the 
choice  and  construction  of  the  numerical  mesh,  that 
is,  the  partitioning  of  the  resolution  domain.  In  the  Fi¬ 
nite  Difference  Method,  this  implies  a  computational 
mesh  with  a  finite  number  of  nodes  where  we  deter¬ 
mine  the  values  of  the  dependent  variables.  Then,  we 
must  approximate  the  existing  derivatives  in  the  gov¬ 
erning  partial  differential  equations. 

Generally,  we  use  two  mesh  systems  in  numerical 
simulation  of  reservoirs:  the  centered  block  system 
and  the  distributed  point  system  [17].  In  this  work,  we 
use  the  former. 

The  mesh  system  contains,  in  the  ^-direction,  rix 
cells  that  we  superimpose  on  the  reservoir.  We  cen¬ 
ter  each  i  cell  on  Xi,  and  we  designate  its  borders  by 
Xi_i/2  and  Xi+i/2.  The  cells  have  dimensions  equal  to 
Axi,  constants  or  not,  satisfying  the  following  relation 

'^Ax,=Ls;.  (35) 

i=l 

SO  that  the  cells  must  cover  the  entire  length  of  the 
reservoir  [1 7].  Similarly,  it  is  also  possible  to  discretize 
the  reservoir  in  the  y-  and  z-  directions. 

The  cells  should  be  small  enough  to  describe 
the  heterogeneous  nature  of  the  reservoir  and,  thus, 
allowing  to  represent  the  flow  characteristics  ade¬ 
quately.  However,  it  is  necessary  to  carefully  deter¬ 
mine  the  total  number  of  cells  in  the  mesh,  because 
the  higher  the  number  of  them,  the  higher  the  num¬ 
ber  of  unknowns  in  the  algebraic  system  that  must  be 
solved,  implying  an  increase  in  computational  effort. 

In  Fig.  1,  a  three-dimensional  reservoir  is  dis¬ 
cretized  into  Ux  blocks  (or  cells)  in  the  x-direction, 
Uy  blocks  in  the  y-direction,  and  blocks  in  the  z- 
direction,  considering  that  nx=ny=nz=3,  for  this  partic¬ 
ular  case.  In  more  general  situations,  the  dimensions 
of  the  blocks  need  not  be  the  same,  and: 

Uy 

^Ayj=Ly,  (36) 

i=i 

riz 

£Azfe=L„  (37) 

k=l 
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where  the  j  and  k  indexes  indicate  the  biocks  in  the 
y-  and  z-  directions,  respectiveiy.  In  this  system,  we 
number  the  border  biocks  by  adding  the  fraction  ±1/2 
to  one  of  the  i,  j,  and  k  indexes,  depending  on  the 
border  that  we  want  to  represent.  For  exampie,  we  in¬ 
dicate  the  boundary  between  ceiis  i,j,  k  and  i  ±  l,j,  k 
by  z±  1/2,/,  k. 


Fig.  1:  Discretized  domain. 


x-direction.  We  can  aiso  obtain  simiiar  forms  for  the  y- 
and  z-  directions,  considering  that  ^  and  Az^  ^  ^ 
are  the  spacing  in  these  respective  directions  and  that 
we  represent  these  boundaries  by  the  indexes  j  ±  1/2 

and  k  ±  1/2. 

Before  proceeding,  we  introduce  transmissibiiities 
as  being  represented  by  T^,,  and  T^: 


yn+l 

■*^^-4-  1  ■  t. 

»±  2  .j,fc 


n+1 


(40) 


.  +  1  ^(AykyY^^ 


(41) 


yn+l 


A,k, 


n+1 


1 

2 


(42) 


and  we  caicuiate  the  properties  at  the  ceii  interface 
using  harmonic  averages  for  rock  and  geometric  prop¬ 
erties,  and  arithmetic  means  for  fluid  properties. 

Again,  empioying  centrai  difference  approxima¬ 
tions. 


n+1 

i+^,j,k 


„"+l 


- 


^i-\-l,j.,k  ^i,j,k 


^i+l,j,fe 

Ax,- 


- 

n%,3,k 


i+l,j,k 


(43) 


3.2.  Numericai  approximation  of  derivatives 

Considering  the  partitioning  iiiustrated  in  Fig.  1 
and  a  totaiiy  impiicit  numericai  scheme,  the  first  term 
containing  the  spatiai  derivatives  of  Eq.  (22)  can  be 
discretized  by  a  centrai  difference  approximation  of 
second-order. 


A 

dx 


n+1 


dx 


i,j,k 


n+1 


i+l/2.j,k 

n+1 


i-ll2,j,k 


(38) 


where 


T,  = 


(39) 


with  simiiar  equations  for  y-  and  z-  directions.  Stiii, 
in  Eq.  (38),  we  do  not  know  the  pressures  at  time 
n  ±  1,  Axij^k  is  the  spacing  of  the  ceii  i  in  the  x- 
direction,  and  the  subscript  i  ±  1/2  indicates  that  the 
variabies  must  be  evaiuated  at  the  ceii  boundary  in  the 


n+1 

i—^,j,k 


^ij,k 


^i—l,j,k 


^i,j.,k  ^i—l.,j,k 


jfi+i 


Ax. 


i—  ^  ,j,k 


(44) 


and  we  can  proceed  simiiariy  to  the  y-  and  z-  direc¬ 
tions. 

Now,  we  appiy  conservative  expansions  to  the 
terms  of  accumuiation  to  preserve  the  mass  baiance. 
It  is  worth  mentioning  that  the  use  of  non-conservative 
schemes  in  the  finite  difference  equations  does  not 
necessariiy  produce  resuits  that  are  not  correct  [17]. 
Therefore,  the  finai  discretized  form  is  given  by: 


4>°C4, 

Bn+1 


i,j,k 


(45) 


where  s,  =  Ax^AyjAzk- 

We  approximate  the  time  derivative  by  a  backward 
Euier  scheme 


(46) 
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Finally,  considering  the  use  of  an  iterative  resolu¬ 
tion  aiming  to  linearize  the  system  of  equations  we 
obtain 


^  ^  •  1  1  ■  u 

/  n+l,-u+l 

\Pi+l,j,k 

n+l,t?+l 

“  P^,j,k 

npn+l,?; 
x .  1  -  . 

f  n+l,v+l 

{Pz-l,j,k 

n+l,t;+l 

“  P^,j,k 

npn+l,?; 

f  n+l,v+l 

\PiJ-\-l,k 

n+l,D+l 

“  P/j,k 

/  n+l,v+l 

{Pi,j-l,k 

n+l.-u+l 

“  Pi,j,k 

ann+l,t; 

f  n+l,v+l 

\PiJ,k-\-l 

n+l,t;+l 
“  P^,J,k 

opn+l,?; 

^  z .  .  ,  1 

/  n+l,v+l 

[Pi,j,k-1 

n+l,T!+l 

“  P^,],k 

) 


) 

) 


-  {Qsc) 


n+l,v+l 

i,j,k 


(Fg) 


n-\-l,v 


+  (^p) 


/  n+l,v+l  _ 

I  Pi,j,k 

At 


(47) 


where  the  level  v  indicates  the  known  values  while  the 
level  -I- 1  the  values  to  be  determined. 

Next,  we  move  on  to  the  discretization  stage  of  the 
well-reservoir  coupling. 


{Qsc) 


n+l,D+l 

id,k 


(  T  n+l,i?+l 


{Pwf) 


n+l,?;+l 

i,j,k 


(48) 


Through  the  term  (Tu,)"^^’’'  the  well-reservoir  cou¬ 
pling  occurs,  following  the  methodologies  based  on 
the  work  of  Peaceman  [26],  which  is  still  the  most 
used  to  date  in  numerical  reservoir  simulation  studies. 

We  note  that  we  use  the  term  to  represent  a 
source  term  for  the  computational  mesh  cell.  In  gen¬ 
eral,  the  well  passes  through  a  set  of  cells  in  the  com¬ 
putational  mesh.  Therefore,  the  expression  to  the  total 
flow  of  the  producing  well  is  given  by: 


Qsc 


Wf 

n+l,v+l 
PiJ,k 

k^Wi 


{Pwf) 


n+l,T!+l 

i,j,k 


(49) 


if  we  assume  a  vertical  well  that  contains  the  cells  Wi 
up  to  Wf  (Fig.  2)  and  we  do  not  consider  the  frictional 
and  convective  effects  inside  the  well. 


Fig.  2:  Vertical  well  in  the  computational  domain. 


We  adopted  as  a  reference  level  k  =  W  tor  calcu¬ 
lating  the  pressure  in  the  well.  Considering  only  the 
gravitational  effects  inside  the  well,  we  can  write: 

/  Nn+1,«+1 

[PwfJijiY 


+ 


where 


E(  J  +  Ti+1, 


k 


k 


(50) 


^i,j,k  —  (P5)ij,fe+l/2  (51) 

and  we  find  the  pressure  in  the  well,  in  each  layer,  in 
the  computational  cell  in  terms  of  wellbore  pressure  in 
the  reference  level  layer: 


(Pwf) 


n+l,'U+l 

i,j,k 


j  Nn+1,D+1 

[PwfJijw 


^payitkhn  -  A.^.W-) ,  (52) 


and  arithmetic  mean  between 

n+l,T! 
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3.3.  Numerical  solution  of  the  algebraic  equation 
system 

The  set  composed  of  Eqs.  (47)  and  (50)  forms 
a  linear  system  of  equations  whose  dependent  vari¬ 
ables  are  the  pressures  in  the  reservoir  and  the  pro¬ 
duction  well. 

We  choose  the  Conjugate  Gradient  Method  [19] 
to  solve  the  resulting  algebraic  equation  system.  Ini¬ 
tially,  it  was  a  direct  method  [19].  However,  it  became 
known  for  its  properties  as  an  iterative  method,  espe¬ 
cially  after  the  development  of  sophisticated  precon¬ 
ditioning  techniques  [32]. 

We  can  increase  the  rate  of  convergence  in  the  it¬ 
erative  process  by  using  preconditioning  techniques 
to  reduce  the  conditioning  number  of  the  matrix  by 
changing  the  original  system  of  equations  and,  thus, 
their  eigenvalues.  Thus,  the  method  used  here  is  the 
Preconditioned  Conjugate  Gradient  (PCG),  using  a  di¬ 
agonal  preconditioning  technique  [7,  17,  32].  If  we 
consider  an  iterative  procedure,  we  attain  the  numeri¬ 
cal  convergence  when 

max  <  tol  (53) 

where  x  represents  the  pressure  in  the  reservoir,  and 
well  and  tol  is  a  numerical  tolerance.  We  use  a  two- 
step  iteration  procedure  [17]: 


4.1 .  Some  well-reservoir  coupling  models 

One  of  the  difficulties  in  modeling  wells  in  field- 
scale  simulations  is  the  fact  that  the  region  where  the 
highest  pressure  gradients  typically  occur  is  close  to 
the  well,  which  is  much  smaller  than  the  usual  di¬ 
mensions  of  a  cell  in  a  computational  mesh  [13].  Ex¬ 
cept  for  simulations  that  somehow  make  use  of  Cylin¬ 
drical  coordinates,  source/sink  terms  are,  in  general, 
used  for  the  implementation  of  internal  boundary  con¬ 
ditions  involving  wells  [3,  17].  Well-reservoir  coupling 
techniques  aim  to  relate  the  pressure  in  the  reservoir 
cells  to  the  pressure  in  the  well,  by  using,  for  exam¬ 
ple,  equations  that  contain  source/sink  terms  and  the 
productivity  index  [35]. 

Among  the  first  works  to  related  the  pressures  in 
the  cell  and  the  well,  we  found  that  of  Van  et  al.  [39]. 
They  stated  that  the  pressure  of  the  cell  containing  the 
well  (for  a  wellbore  centered  in  it)  must  be  equal  to  the 
average  volumetric  pressure  of  the  cell.  Thus,  for  the 
hypotheses  of  steady-state  flow,  the  average  pressure 
p  (or  the  pressure  in  the  cell),  located  in  the  reservoir 
between  the  radii  and  ri,  (n  »  I'w)  is  given  by: 


where  n  is  called  the  block  radius  and  determined,  for 
a  vertical  well,  as 


P  =  Pwf  + 


2'nqscBpL 

kh 


1 .  in  a  set  of  internal  iterations,  we  apply  the  PCG 
method  to  obtain  the  pressures,  and 

2.  in  a  set  of  external  iterations,  we  update  the  co¬ 
efficients,  such  as  the  transmissibilities. 

IV.  WELL-RESERVOIR  COUPLING 

The  modeling  of  wells  in  the  numerical  simula¬ 
tion  of  reservoirs  presents  some  difficulties  as,  for 
example,  the  discrepancies  in  spatial  scales  when 
comparing  the  wellbore  radius  with  the  dimensions 
of  the  oil  reservoir,  the  geometry  of  complex  wells, 
the  flow  transition  from  the  porous  medium  to  a  duct, 
among  other  factors.  In  this  section,  we  investigate 
the  treatment  of  the  term  source,  review  some  well- 
reservoir  coupling  models,  focusing  on  the  Peace- 
man  [26]  model  and  its  extensions,  and  discuss  the 
methodology  that  incorporates  the  transient  effects, 
both  in  the  productivity  index  and  in  the  equivalent  ra¬ 
dius. 


We  observe,  then,  that  a  methodology  to  deal  with 
the  problem  of  the  representation  of  wells  in  Cartesian 
coordinates  consists  of  initially  considering  the  one¬ 
dimensional  and  radial  mass  flow  around  the  well  [3]. 
The  idea  is  to  reconcile  the  analytical  solutions,  which 
we  can  obtain  for  the  one-dimensional  flow  in  Cylindri¬ 
cal  coordinates,  with  the  discretized  equations,  written 
in  Cartesian  coordinates. 

4.2.  Peaceman’s  equivalent  radius 

We  can  obtain  the  first  models  for  well-reservoir 
coupling  by  considering  a  two-dimensional  flow  gov¬ 
erned  by  the  Darcy  equation  written  in  Cartesian  co¬ 
ordinates  and  a  one-dimensional  flow  in  a  cylindrical 
reservoir.  Peaceman  [26],  to  obtain  its  well-reservoir 
coupling  model,  used  the  discrete  equation  for  the 
pressure  p{x,y)  in  the  two-dimensional  flow  together 
with  the  analytical  solution  for  the  pressure  p{r)  in  the 
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one-dimensional  flow,  both  for  a  porous  medium  of  the 
same  length  Az. 

From  the  hypotheses  of  steady-state  flow,  homo¬ 
geneous  and  isotropic  porous  medium  (kj,  -  ky  ^ 
kn),  assuming  that  Ax  =  Ay,  and  symmetry,  we  ob¬ 
tain  from  the  discrete  form  of  HDE  [17] 


Pi+l,j  —  Pi,j 


Byqsc 

AknAz 


(56) 


On  the  other  hand,  we  know  that  the  analytical 
solution  of  the  one-dimensional  flow  in  the  radial  di¬ 
rection,  in  a  reservoir  described  in  Cylindrical  coordi¬ 
nates,  for  a  steady-state  flow,  is  given  by  [12] 


P  (r)  =  Pwf  - 


qscBp 
77—. — ^  in 
2'KkH^Z 


(57) 


Peaceman  [26]  considered  the  pressure  p{r)  coin¬ 
ciding  with  the  pressure  Pi+i,3-  So,  for  r  =  Ax,  we 
obtain 


p  (Ax)  =  pi+ij  =  pi,f  - 


qscBp  . 

„  ,  ^  In 
27rfcAz 


(58) 


Combining  Eqs.  (56)  and  (58),  that  is,  consider¬ 
ing  the  discrete  form  in  Cartesian  coordinates  and 
the  analytical  expression  for  the  Cylindrical  geometry 
around  the  well,  we  can  arrive  in  an  equation  for  the 
flow  in  the  well,  written  using  the  productivity  index. 


Qsc 


27rfcgAz(p^j  -Pwf) 


Gw  , 
Jw  {jPi,j 


Pwf) 
Pwf)  j 


(59) 


where  re,  is  Peaceman’s  equivalent  radius,  re,  « 
0.198Aa;  [17], 

Since  then,  researchers  have  created  several 
extensions  assuming  a  steady-state  flow,  from 
the  Peaceman  model  for  the  equivalent  radius  [26], 
to  improve  the  numerical  approximation  considering 
other  hypotheses.  For  example,  the  concept  of  the 
equivalent  radius  was  investigated  for  the  case  of 
fca;  ^  ky  and  Ax  ^  Ay,  thus  allowing  us  to  address 
more  comprehensive  situations  in  numerical  simula¬ 
tions  [27],  Using  coordinate  transformations.  Peace- 
man  [27]  obtained  the  following  expression  for  the 
equivalent  radius 


which  we  must  use  in  conjunction  with  the  equation 


2,7:  s/JT/Ty  A  z 

= - /„  N  (Pi,j  -  Pwf)  ■ 


(61) 


It  is  worth  mentioning  that  this  derivation  for  the  equiv¬ 
alent  radius,  for  an  anisotropic  medium,  considers  that 
Ax  ^  Ay  but  that  the  mesh  is  uniform  and  that  the 
well  is  far  from  the  reservoir  borders. 

Other  extensions  have  been  developed,  such  as 
the  one  used  by  Al-Mohannadi  et  al.  [2] 


'Tj/;  — 


2Try/kykzAx 


Ax 


Bp 


'-(%)1>"(S) 


(62) 


where 


(63) 


for  a  horizontal  well  parallel  to  the  x-axis. 

Another  well-coupling  model,  also  for  a  horizontal 
well,  is  given  by  [4,  5] 


2'K^JkykzAx 


where 


Teq  =  Q.l^ikykx) 


1/4  ( A^  Az^ 

ky  ky 


2\  1/2 


1  +  exp 


2.215-3 


QC  f  '^y'^y  A 

V  / 


1  -F  0.533  ^ 


(¥) 


(64) 


(65) 


for  a  well  parallel  to  the  x-axis,  where  an  = 
{Ay/ Ax)  y^kz/ky.  As  examples  of  other  methodolo¬ 
gies  we  can  quote  Peaceman  [28],  Peaceman  [29], 
and  Babu  et  al.  [6].  Researchers  have  even  studied 
inclined  wells  [22]. 

So  that  we  can  apply  these  well-reservoir  cou¬ 
pling  models,  some  restrictions  must  be  considered, 
such  as,  for  example,  the  existence  of  a  minimum  dis¬ 
tance  between  the  wells  and  the  reservoir  borders, 
and  between  the  wells  themselves.  Despite  its  limi¬ 
tations  [23],  we  extensively  use  the  technique  intro¬ 
duced  by  Peaceman  [26]  in  the  numerical  simulation 
of  reservoirs. 
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4.3.  Extensions  for  the  transient  productivity  index 
As  aiready  said,  the  resuits  obtained  with  the  tech¬ 
nique  suggested  by  Peaceman  [26]  are  subject  to  a 
numericai  artifact,  which  is  not  reiated  to  the  phe¬ 
nomenon  of  physicai  weiibore  storage.  Fundamen- 
taiiy,  this  happens  due  to  the  use  of  the  hypothe¬ 
sis  of  steady-state  flow  near  the  weii.  In  reality,  we 
should  have  considered  the  fluid  flow  as  being  tran¬ 
sient.  Therefore,  other  models  have  been  developed, 
including  the  well-reservoir  coupling  for  the  transient 
regime.  For  flow  in  the  transient  and  pseudo-steady 
state  regime,  in  the  vicinity  of  the  well,  Peaceman  [26] 
suggested  that  we  can  obtain  the  equivalent  radius 
from  the  following  equation: 

Teq  =  [4f_D  exp(-A  -  4ttpd)]AL  (66) 


Peaceman  technique  for  the  initial  times,  but  it  can 
present  deviations  when  the  transient  regime  appears 
in  the  porous  medium  [9]. 

Therefore,  instead  of  using  the  equivalent  radius 
as  req  =  0.198  Aa;,  it  is  possible  to  calculate  its  tran¬ 
sient  version.  The  pressure  variation,  for  radial  flow  in 
the  transient  regime  and  a  vertical  well,  is  given  by  [2]: 


Ap  =  — 


Ln^/kxky  ^  \ty/k^k. 


qscBp  ^  I  (j^ctprl^ 

\  +  /  ly> 

y  > 


(70) 


From  Eq.  (70),  we  can  use  the  Newton-Raphson 
method  to  determine  the  req  as  a  function  of  time. 


r’e(j„+i  —  req^ 


fjreqj 

f'il-eqS 


(71) 


where 

is  the  dimensionless  time  and  AL  {AL  =  Ax  =  Ay)  is 
the  spatial  increment  of  the  computational  mesh,  and 


PD  =  -  (  )  {Pini  -  Pt,3,k)  (68) 

\qscBfiJ 

is  the  dimensionless  pressure,  and  p,^j^k  is  the  pres¬ 
sure  in  the  cell  that  contains  the  well.  Blanc  et  al.  [9] 
also  concluded  that  when  to  >  '\,  the  equivalent  ra¬ 
dius  calculated  using  Eq.  (66)  is  close  enough  to  its 
value  determined  by  Peaceman’s  model  for  steady- 
state  flow. 

Over  time,  we  have  been  using  analytical  solu¬ 
tions  capturing  the  transient  regime,  in  the  porous 
medium,  to  interpret  the  results  of  pressure  tests  in 
wells.  Theis  [37],  for  example,  offered  one  of  the  first 
solutions  to  the  problem  of  pressure  drop  in  a  well. 

Blanc  et  al.  [9]  presented  a  correction  to  the  orig¬ 
inal  Peaceman  model  and  introduced  the  Transient 
Well  Index: 


Jtwi 


pB 


El 


Airkh 

(  (t>ctprl\  ^  ((t>pctr%\ 


where  the  subscripts  n  and  n  +  1  represent  the  itera¬ 
tions  in  the  Newton-Raphson  method  and 


fireqj  =  - 


qscBp 


AnAzyJkxk., 


-.El 


eqn 


\Jkx  ky 

(fipct 


-  Ap, 


(72) 


while  its  derivative  is  given  by 


fireqj 


qscBp 


2TTreq„Az,Jkxk, 


•  exp 


(  \/  kx  ky 

y  (t>pct 


t 

(73) 


In  Eq.  (72),  Ap  =  p™  where  p,™  is  the  initial 

pressure  of  the  reservoir,  and  p,^j^k  is  the  pressure  in 
the  cell  containing  the  well. 

From  these  equations,  we  can  obtain  the  expres¬ 
sion  for  the  determination  of  the  equivalent  transient 
radius 


qscBpEi(a)  AApnAz^Jkxky 
2qBp  exp(— cr) 


(74) 


(69)  where 


where  Ei{u)  =  -Ei{-u)  and  req  =  0.198  AL  (steady- 
state  flow  hypothesis).  This  model  offers  better  so¬ 
lutions  than  those  obtained  with  the  conventional 


tj  = 


-s/kxky 

4>CtP 


(75) 
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and 


—  Pini  Pi,j,k- 


(76) 


In  this  work,  p,j^k  is  evaluated  at  time 

When  the  boundary  effects  start  to  affect  the  well¬ 
bore  pressure  behavior,  re,(f)  becomes  too  large  to 
represent  the  flow  dynamics.  In  this  case,  we  calcu¬ 
late  Teqit)  from  Eqs.  (70),  (71)  and  (72)  while  t  <  t*, 
where  t*  is  the  characteristic  time  associated  with 
the  boundary  effects  [2],  For  t  >  t*,  we  must  use 
req  =  regit*),  and  we  keep  the  transient  equivalent  ra¬ 
dius,  calculated  at  t*,  frozen  until  the  end  of  the  simu¬ 
lations. 

Here,  a  criterion  to  shift  the  equivalent  transient  ra¬ 
dius  calculation  is  considered,  based  on  the  time  re¬ 
quired  to  reach  the  pseudo-steady  regime.  Therefore, 
we  use  the  concept  of  dimensionless  time  considering 
A  (the  area  normal  to  the  z-axis)  as  a  reference  [30] 
so  that 


kt 

(jipLCtA  ’ 


(77) 


and  we  assume  that  the  producing  well  is  in  the  center 
of  the  reservoir  and  that  it  has  a  square  drainage  area. 
For  toA  <  0.09,  the  reservoir  still  behaves  as  being  in¬ 
finite  (transient  regime).  On  the  other  hand,  for  tDA= 
0.1 ,  we  reach  the  pseudo-steady  flow  regime  [30]. 


V.  NUMERICAL  RESULTS 

In  the  study  carried  out  in  this  work,  the  well- 
reservoir  coupling  models,  including  transient  effects, 
were  incorporated  into  our  simulator,  developed  in  C 
programming  language. 

We  calculate  the  well  pressure  values  considering 
a  constant  production  flow  rate  under  standard  condi¬ 
tions.  We  use  two  types  of  graphs  in  the  analysis  of 
the  results: 

1 .  specialized  plot:  well  pressure  curve  as  a  func¬ 
tion  of  time; 

2.  diagnostic  plot:  pressure  drop  curves  in  the  well, 
Apu,/,  and  the  Bourdet  derivative  [11] 

as  a  function  of  time. 

Our  approach  is  in  the  context  of  reservoir  simula¬ 
tion  and  the  transient  analysis  of  pressure  tests.  We 


know  that  from  the  pressure  values,  it  is  possible  to 
determine  some  reservoir  properties,  useful  for  plan¬ 
ning  the  completion  of  the  well  or  the  reservoir  de¬ 
pletion  plan.  In  addition  to  the  well-reservoir  coupling 
model  of  Peaceman  [26],  considering  a  steady-state 
flow,  three  extensions  were  implemented,  assuming  a 
transient  flow  regime.  In  all  simulations,  we  stipulate 
that  Ax  =  Ay  and  =  ky  =  =  k. 


5.1 .  Model  1 

It  is  the  conventional  model  of  Peaceman  [26], 
where  we  assume  a  steady-state  flow  and  whose  pro¬ 
ductivity  index  and  the  equivalent  radius  considered 
are,  respectively,  given  by 


2'KkAz 


(79) 


and 


Teq  =  0.198Ax. 


(80) 


5.2.  Model  2 

Peaceman  [26]  also  proposed  this  model.  How¬ 
ever,  it  considers  that  the  flow  is  transient  in  the  cells 
containing  the  well.  The  productivity  index  is  the  same 
as  for  Model  1,  although  the  equivalent  radius  incor¬ 
porates  the  transient  effects: 


2'KkAz 


(81) 


where 

req  =  [4t_D  exp(-A  -  4kpd)]^^^Ax,  (82) 

A  =  exp(0. 57722),  0.57722  is  Fuler’s  constant, 

/  \ 

PD  =  -  (  - ^  {Pxni  -  Pi,j,k)  ,  (84) 

and  we  use  here  the  criterion  provided,  based  on  toA, 
to  shift  the  calculation  of  the  equivalent  radius  (we  do 
not  update  r^q  when  border  effects  occur). 
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5.3.  Model  3 

Blanc  et  al.  [9]  originally  proposed  this  model, 
and  it  employs  a  transient  productivity  index.  Nev¬ 
ertheless,  the  equivalent  radius  proposed  by  Peace- 
man  [26],  for  the  steady-state  flow  (Model  1),  is  ap¬ 
plied 


and  high  viscosities,  using  power  series  to  approxi¬ 
mate  the  exponential  integrals  is  not  suitable  for 
us.  Therefore,  we  chose  the  approximation  of  the  ex¬ 
ponential  integral  Ei{x)  suggested  by  Segletes  [34], 
based  on  a  polynomial  fit  [31].  It  is  not  our  knowledge 
that  others  have  done  this  previously. 


Jtwi 


47rfcAz 


p,B 


El  -  El 


/  (pctiirl^  \ 
4kt  ) 


where 


(85) 


Teg  =  O.lQSAx.  (86) 

5.4.  Model  4 

In  this  well-reservoir  coupling  model,  we  consider 
the  incorporation  of  transient  effects  in  the  productivity 
index  and  the  equivalent  radius  [2,  9].  The  productiv¬ 
ity  index  is  the  same  as  that  defined  in  Model  3,  and 
we  calculate  the  equivalent  transient  radius  using  the 
Newton-Raphson  method.  Thus,  we  have: 


Jtwi  = - 

p.B 

knowing  that 

^69.14 

and 


AirkAz 
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fjreqj 

fireqj" 


\  4fci 
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—  Ap. 


(87) 
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(89) 


The  criterion  based  on  Eq.  (77)  is  also  applied  to  in¬ 
terrupt  the  update  of  the  transient  equivalent  radius. 

Many  analytical  approaches  to  evaluate  Ei{u) 
have  appeared  in  the  literature.  Tseng  and  Lee  [38] 
provide  an  excellent  survey  on  the  subject.  Note  that 
there  is  no  single  analytical  approximation  that  is  valid 
for  the  entire  range  for  u  >  0.  Tseng  and  Lee  [38] 
claim  that  different  methods  are  suitable  for  approxi¬ 
mating  Ei{u)  for  a  wide  range  of  u  values.  Neverthe¬ 
less,  in  practical  applications,  we  often  do  not  have 
only  a  single  calculation  method  that  covers  the  entire 
range  of  interest  [8]. 

In  this  work,  as  we  are  interested  in  the  short¬ 
est  times  and  flows  that  may  have  low  permeabilities 


5.5.  Analytical  solution  and  mesh  refinement 

Now,  we  compare  the  results  obtained  with  the 
different  models  of  well-reservoir  coupling.  Table  1 
contains  the  parameters  used  in  the  construction  of 
the  standard  simulation  case.  Unless  when  explicitly 
mentioned,  the  data  in  this  table  are  those  used  in  all 
simulations  performed. 

Table  1 :  Parameters  for  standard  simulation. 


Parameter 

Value 

Unit 

B° 

1.3 

RB/STB 

Co 

4x10-® 

psi-^ 

Cf, 

2x10-® 

psi-^ 

3x10-® 

psi-^ 

kx,  ky  and  kz 

1x10-2 

Darcy 

Lx 

10,000 

ft 

Ly 

10,000 

ft 

Lz 

80 

ft 

Lwf 

80 

ft 

rix 

81 

- 

Uy 

81 

- 

riz 

5 

- 

Pint 

8,000 

psi 

pO 

8,000 

psi 

Q  sc 

-400 

STB/day 

Tw 

0.2 

ft 

tol 

1x10-® 

psi 

^max 

730 

day 

5  At 

1.2 

- 

^^ini 

1x10-® 

day 

^^max 

10 

day 

A 

1.1 

cp 

P 

52.4 

Ib/ft  ® 

(t> 

0.2 

- 

(j)^ 

0.2 

- 

In  this  table,  L^,/  is  the  length  of  the  producing 
well,  tmax  represents  the  maximum  production  time, 
6At  is  the  factor  used  to  vary  the  time  increment 
(At"+^  =  and  Atin,  and  Atmax  are  the  initial 
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and  maximum  values  of  the  time  increment,  respec¬ 
tively. 

Here,  the  results  show  the  pressure  curves  in 
the  well,  the  pressure  drop,  and  the  Bourdet  deriva¬ 
tive  [11,  10]  as  a  function  of  the  elapsed  production 
time. 

The  first  test  performed  was  done  to  compare  re¬ 
sults  from  well-reservoir  coupling  models  with  an  an¬ 
alytical  solution  for  the  particular  case  of  constant  vis¬ 
cosity  and  formation  volume  factor.  Figures  3  and  4 
show  comparisons  of  the  results  obtained  with  Mod¬ 
els  1  and  4,  using  the  standard  set  of  parameters  for 
the  simulations,  with  the  results  obtained  with  the  an¬ 
alytical  solution  given  by  Ozkan  [25]. 

It  should  be  made  clear  that  the  analytical  solu¬ 
tion,  as  already  said,  makes  use  of  simplifying  as¬ 
sumptions,  such  as  B  and  fj.  constants.  Besides,  we 
must  apply  a  numerical  procedure  to  obtain  the  ana¬ 
lytical  solution,  which  depends  on  the  evaluation  of 
functions  of  Bessel  [24]  and  numerical  inversion  of 
Laplace  transform  [25],  performed  using  the  Stehfest 
algorithm  [36].  Therefore,  the  results  for  the  analytical 
solution  are  not  available  in  the  entire  range  of  nu¬ 
merical  results.  Even  so,  it  is  possible  to  observe  that 
the  results  are  in  good  agreement  with  the  numerical 
results  obtained  using  Model  4.  On  the  other  hand, 
we  note  a  deviation  when  there  is  a  comparison  with 
Model  1  (Fig.  3).  When  the  border  effects  are  present, 
the  three  solutions  have  very  similar  behaviors  (about 
250  days  of  production).  Figure  4  corroborates  the 
discussions  about  Fig.  3,  and  we  see  that  the  Bour¬ 
det  derivative  highlights  the  difference  in  the  behavior 
of  the  numerical  solution  and  the  effect  of  the  numer¬ 
ical  artifact  due  to  the  use  of  a  constant  equivalent 
radius  (r^,). 

Also  noteworthy  is  the  qualitative  similarity  of  the 
result  obtained  with  the  original  technique  of  Peace- 
man  [26]  (Model  1)  with  the  behavior  of  the  pres¬ 
sure  when  there  is  physical  storage  in  the  well.  Also, 
there  is  a  discontinuity  in  the  derivative  curves  corre¬ 
sponding  to  the  transient  and  the  pseudo-permanent 
regimes.  It  is  due  to  the  strategy  employed  to  change 
the  calculation  of  the  equivalent  radius.  From  these 
results,  it  is  possible  to  conclude  that  the  results  of 
Model  4  present  a  behavior  compatible  with  that  ex¬ 
pected  for  the  real  problem. 

Aiming  to  compare  the  four  well-reservoir  coupling 
models  that  we  implemented  in  this  work,  in  Fig.  5, 


we  show  the  set  of  results  for  the  pressure  in  the  well. 
We  can  note  the  existence  of  a  plateau  in  the  curve 
obtained  with  Model  1  for  the  initial  production  times. 
For  the  standard  case,  this  plateau  (numerical  stor¬ 
age)  lasts  approximately  one  day  of  production,  inter¬ 
fering  in  the  pressure  response,  reaching  a  difference 
of  about  33  psi  for  the  same  time  when  compared  to 
the  Model  4  (in  principle,  the  most  correct).  Regarding 
the  models  that  incorporate  the  transient  effects  in  the 
equivalent  radius.  Models  2  [26]  and  4,  we  observe 
that  we  obtained  better  results  with  these  models  (su¬ 
perimposed  in  Fig.  5)  in  comparison  with  those  based 
on  the  steady-state  flow  assumption.  Besides  that, 
in  the  initial  instants  of  production,  the  numerical  arti¬ 
fact  no  longer  appears,  and  we  correctly  capture  the 
expected  flow  regimes.  However,  according  to  Blanc 
et  al.  [9],  the  results  of  Model  2  are  not  always  as  fa¬ 
vorable  as  those  obtained  here.  For  example,  it  can 
happen  if  we  vary  the  reservoir  and  fluid  properties. 

Given  the  results  obtained  with  Model  3  (Fig.  5), 
it  is  possible  to  observe  the  numerical  artifact  occur¬ 
ring.  However,  with  less  intensity  than  that  of  Model  1 , 
and  we  do  not  perceive  it  when  the  production  begins. 
It  appears  approximately  15  minutes  after  the  begin¬ 
ning  of  the  production.  Nevertheless,  after  about  one 
day,  the  results  show  behavior  consistent  with  those 
of  Models  2  and  4.  We  should  remark  that  the  re¬ 
sults  of  Model  3  are  in  line  with  that  reported  by  Blanc 
et  al.  [9].  On  the  other  hand,  in  Model  4,  when  we 
incorporate  transient  effects  in  both  productivity  index 
and  equivalent  radius,  the  numerical  artifact  does  not 
exist  anymore,  and  the  results  are  physically  correct 
for  the  entire  production  duration.  For  the  four  mod¬ 
els,  we  capture  the  border  effect  when  we  reach  a 
time  of  approximately  250  days  of  production.  Among 
the  four  models,  we  can  say  that  Models  2  and  4  were 
the  ones  that  presented  the  best  results  considering 
the  transient  flow  regime.  Finally,  we  should  stress 
that  the  behavior  of  Model  4  results  is  also  in  line  with 
those  reported  by  Blanc  et  al.  [9]. 

Furthermore,  we  also  carried  out  a  mesh  refine¬ 
ment  study  using  Models  1  and  4  that  take  into  ac¬ 
count  the  transient  effects  of  flow.  Table  2  shows  the 
number  of  cells  that  we  used  in  the  generation  of  the 
different  meshes  that  we  employed  in  the  mesh  refine¬ 
ment  study.  Here,  Ux,  Uy,  and  are  the  number  of 
cells  used  to  discretize  the  oil  reservoir  in  the  x-,  y- 
and  z-  directions,  respectively. 
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4 


Fig .  3 :  Comparison  of  the  results  of  Models  1  and  4  with  the  analytical  solution.  Specialized  plot. 


4 


Fig.  4:  Comparison  of  the  results  of  Models  1  and  4  with  the  analytical  solution.  Diagnostic  plot. 
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Fig.  5:  Comparison  of  the  results  of  the  four  well-reservoir  coupling  models. 


We  can  see  the  results,  well  pressure  as  a  function 
of  time,  in  Figs.  6  and  7. 


Table  2:  Meshes. 


Mesh 

rix 

Uy 

n-z 

1 

11 

11 

5 

2 

21 

21 

5 

3 

41 

41 

5 

4 

81 

81 

5 

5 

161 

161 

5 

For  Model  1 ,  it  was  possible  to  observe  a  numeri¬ 
cal  convergence  as  we  refine  the  mesh.  Except  for  the 
initial  moments,  where  the  solutions  are  quite  differ¬ 
ent.  We  can  also  realize  that  the  mesh  refinement  in¬ 
fluences  the  magnitude  of  the  numerical  artifact,  since 
the  more  refined  the  mesh,  the  smaller  the  storage 
phenomenon.  For  example,  for  Mesh  1  {nx=ny='\'\), 
we  noted  that  the  numerical  artifact  lasted  approxi¬ 
mately  30  days,  while  for  Mesh  5  (nx=ny='\Q'\)  the 
duration  was  about  0.1  days. 

On  the  other  hand,  we  did  not  observe  the  same 
behavior  when  we  used  Model  4  (Fig.  7),  which  incor¬ 


porates  the  transient  effects.  In  this  case,  the  pres¬ 
sure  variation  presents  a  behavior  compatible  with 
that  of  real  fluid  flow  in  an  oil  reservoir.  Although  the 
results  are  practically  overlapping,  for  all  meshes  and 
time,  including  the  border  effects,  numerical  conver¬ 
gence  did  not  occur  in  a  similar  way  to  that  of  Model  1 . 
However,  the  difference  between  the  values  is  small, 
and  we  cannot  see  it  in  the  graph.  It  is  worth  point¬ 
ing  out  that  the  same  issue  was  reported  in  the  lit¬ 
erature  [1,  2],  indicating  a  possible  loss  of  accuracy 
when  we  utilize  very  refined  meshes.  However,  this 
is  still  an  open  problem.  Another  important  conclu¬ 
sion  that  we  can  draw  is  the  fact  that  we  can  use  less 
refined  meshes,  because  the  difference  between  the 
pressure  values,  obtained  with  more  refined  meshes, 
is  almost  imperceptible.  Therefore,  this  fact  implies 
less  computational  effort. 

We  also  analyzed  the  influence,  on  well  pressure 
variation,  of  the  growth  rate  of  the  time  step  in  the 
well-reservoir  coupling  for  Model  4.  We  obtained  the 
results  considering  three  different  values  of  the  growth 
rate  of  the  time  step,  (5At=1-15,  1.20,  and  1.25.  We 
show  the  results  in  Figs.  8  and  9,  for  the  pressure  in 
the  well,  the  pressure  drop,  and  the  derivative  of  the 
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Fig.  6:  Result  of  mesh  refinement  for  Mode!  1. 
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Fig.  7:  Result  of  mesh  refinement  for  Model  4. 
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pressure  drop.  For  all  the  values,  we  noticed  that  the 
results  are  practically  the  same.  Furthermore,  they 
coincide  throughout  all  simulation  time  and  present 
the  behavior  expected  for  this  type  of  flow.  That  Is,  for 
the  single-phase  flow  of  oil  In  a  reservoir  containing  a 
vertical  producer  well. 

5.6.  Low  permeability  and  high  viscosity 

We  also  performed  other  simulations  varying  the 
permeability  and  viscosity.  However,  we  employed 
lower  values  for  permeability  and  higher  for  viscos¬ 
ity.  In  this  study,  we  only  considered  Models  1  and  4 
and  specialized  plots.  We  did  a  more  detailed  analy¬ 
sis  to  compare  the  two  well-reservoir  coupling  models, 
for  these  permeability  and  viscosity  values,  to  under¬ 
stand  how  the  choice  of  the  model  may  Impact  the 
results  when  the  numerical  artifact  Is  more  significant. 
We  considered  a  total  production  time  of  80  days  and 
a  prescribed  flow  (qsc)  of  -100  STB/day.  This  ap¬ 
proach  alms  to  deal  with  situations  closer  to  those  of 
unconventional  reservoirs.  This  type  of  study  has  be¬ 
come  a  trend  more  recently,  and.  In  general,  the  effect 
of  numerical  storage  tends  to  be  more  pronounced. 

For  the  model  of  Peaceman  [26]  (Model  1),  It  Is 
possible  to  observe  that  the  permeability  value  di¬ 
rectly  Influences  the  magnitude  of  the  numerical  ar¬ 
tifact  In  the  Initial  moments.  The  higher  the  perme¬ 
ability,  the  smaller  the  size  of  the  artifact  (Fig.  10).  It  Is 
physically  consistent  because  when  the  permeability 
Is  higher,  the  lower  the  flow  resistance  In  the  porous 
medium.  Therefore,  It  reaches  the  transient  regime  In 
the  porous  medium  more  quickly,  as  well  as  for  border 
effects.  It  Is  worth  noting  that  since  permeability  has  a 
value  considered  low  for  this  model  of  a  well-reservoir 
system,  the  numerical  artifact  has  a  longer  duration. 
Thus,  If  we  Intend  to  study  the  phenomena  that  occur 
In  the  Initial  moments  (or  even  for  longer  times)  with 
this  model,  the  results  will  not  be  accurate. 

As  an  Illustration,  In  Fig.  10,  for  the  lowest  per¬ 
meability  (fc=0.25x10“^  Darcy),  the  numerical  arti¬ 
fact  Impacted  almost  all  results.  The  pressure  differ¬ 
ence,  comparing  the  two  models,  was  376  psi  at  t=2 
days  of  production  and  54.6  psI  at  t=70  days.  For 
A:=1. 00x10'^  Darcy,  the  magnitude  of  the  numerical 
artifact  was  smaller,  but  not  negligible  since,  after  0.6 
and  19  days  of  production,  there  was  a  pressure  dif¬ 
ference  of  about  94  and  1 0  psi,  respectively. 

It  was  also  possible  to  observe  that  for  these  low 


values  of  permeability,  besides  the  duration  of  the  nu¬ 
merical  artifact  being  longer.  In  some  cases.  It  was 
not  possible  to  capture  the  border  effects.  It  Is  rele¬ 
vant  In  the  context  of  the  analysis  of  pressure  tests, 
as  It  Indicates  the  real  need  for  long-duration  tests 
when  It  Is  necessary  to  determine  boundary  effects. 
We  know  that  It  Is  due  to  the  great  difficulty  of  the  fluid 
to  flow  through  the  porous  medium  and,  therefore.  It 
will  take  longer  to  perceive  the  reservoir  boundary  ef¬ 
fects  (Fig.  10).  Therefore,  we  can  conclude  that  the 
transient  well-reservoir  coupling  model  presented  the 
best  results  since  we  were  able  to  prevent  the  appear¬ 
ance  of  the  numerical  artifact.  Besides  that,  the  use  of 
the  expression  of  Segletes  [34]  to  calculate  Ei  allows 
us  to  obtain  results  In  a  range  In  which  other  formulas 
failed. 

Contrary  to  what  happens  for  the  permeability  vari¬ 
ation,  for  the  conventional  model  [26]  (Model  1),  the 
numerical  artifact  has  a  longer  duration  and  magni¬ 
tude  as  the  viscosity  Increases  (Fig.  11).  Therefore, 
the  time  of  occurrence  and  the  magnitude  of  the  ar¬ 
tifact  vary  depending  on  the  viscosity.  Outside  the 
region  associated  with  the  numerical  artifact,  the  be¬ 
havior  of  the  results  of  the  conventional  model  quali¬ 
tatively  tends  to  be  physically  correct.  Also,  we  know 
that  border  effects  occur  more  quickly  for  lower  viscos¬ 
ity  values.  Regardless,  we  did  not  detect  the  bound¬ 
ary  effects  for  any  of  the  viscosity  values  proposed  In 
our  test. 

VI.  CONCLUSION 

The  objective  of  this  work  was  to  Implement  a 
model  for  the  transient  well-reservoir  coupling  In  a 
reservoir  simulator  to  correctly  describe  the  pressure 
behavior  In  the  well  In  the  Initial  production  times.  We 
also  Investigated  the  effects  of  the  mesh  refinement 
and  the  size  of  the  time  step  on  the  results. 

As  well  known,  the  traditional  well-reservoir  cou¬ 
pling  model  of  Peaceman  [26]  Is  still  widely  used 
nowadays.  However,  we  saw  that  It  has  a  flaw  that 
leads  to  numerical  storage  at  the  beginning  of  the 
simulation.  Nevertheless,  It  Is  simple,  and  we  can  cor¬ 
rectly  determine  the  wellbore  pressure  In  a  wide  range 
of  applications  except  for  the  Initial  Instants  of  produc¬ 
tion.  We  must  also  emphasize  that  we  were  able  to 
detect  the  border  effects  with  this  technique,  without 
the  need  for  a  specific  criterion  to  change  the  calcula¬ 
tion  of  the  equivalent  radius. 
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Fig.  8:  Results  for  different  S At  growth  ratios.  Specialized  plot. 
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Fig.  9:  Results  for  different  6At  growth  ratios.  Diagnostic  plot. 
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Fig.  10:  Results  for  different  low  permeabilities.  Specialized  plot 
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Fig.  11:  Results  for  different  high  viscosities.  Specialized  plot. 
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In  addition  to  physicai  properties,  the  numericai  ar¬ 
tifact  is  aiso  affected  by  mesh  refinement.  The  refine¬ 
ment  can  be  used,  in  some  cases,  to  mitigate  the  non- 
physicai  storage,  but  it  ieads  to  an  increase  in  compu- 
tationai  cost.  On  the  other  hand,  despite  the  quaii- 
tativeiy  and  quantitativeiy  correct  resuits  that  we  ob¬ 
tained  with  Modei  4,  its  numericai  convergence  must 
be  better  understood,  as  aiready  pointed  out  by  other 
authors.  However,  this  is  not  a  big  probiem,  because 
we  have  not  to  use  refined  computationai  meshes 
when  considering  Modei  4,  as  we  couid  see  in  this 
work.  Nevertheiess,  this  issue  deserves  further  stud¬ 
ies. 

We  must  aiso  stress  that  the  use  of  the  expression 
proposed  by  Segietes  [34]  to  caicuiate  the  exponentiai 
integrai,  aiiowed  us  to  consider  flows  with  iow  perme- 
abiiities,  high  viscosities,  and  short  production  time. 
Sometimes,  this  is  not  possibie  when  we  appiy  other 
expressions  avaiiabie  in  the  iiterature. 
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